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Abstract 

This paper is a continuation of our earlier work [3] in which a numerical moment 
method with arbitrary order of moments was presented. However, the computation 
may break down during the calculation of the structure of a shock wave with Mach 
number Mq ^3. In this paper, we concentrate on the regularization of the moment 
systems. First, we apply the Maxwell iteration to the infinite moment system and 
determine the magnitude of each moment with respect to the Knudsen number. After 
that, we obtain the approximation of high order moments and close the moment sys- 
tems by dropping some high-order terms. Linearization is then performed to obtain 
a very simple regularization term, thus it is very convenient for numerical implemen- 
tation. To validate the new regularization, the shock structures of low order systems 
are computed with different shock Mach numbers. 

Keywords: Boltzmann-BGK equation; Maxwellian iteration; Regularized moment 
equations 

1 Introduction 

In the field such as high altitude flight and microscopic flows, gas is considered to be 
very rarefied and outside the hydrodynamic regime. In this case, usual fluid models such 
as Euler equations and Navier-Stokes-Fourier system will fail when the rarefied effect is 
significant. The moment method, which was first proposed by Grad [6], is focused on 
the description of the rarefied gases using a small number of variables. Almost all mo- 
ment methods are derived from the Boltzmann equation which is regarded to be able to 
capture the rarefied effects accurately. In [3] , a special expansion of the distribution func- 
tions is adopted to make it possible the solve the associated Grad-type moment equations 
numerically without the explicit expressions of the system. Then, we followed [16] and 
numerically regularized the system using the technique of a modified Chapman-Enskog 
expansion. In [3], it has been verified numerically that a smooth shock structure with 
Mach number M$ = 2 can be obtained by solving the R20 equations with a Riemann 
problem until the steady state. However, it was found in our numerical experiments that 
if we set the shock Mach number Mq ^ 3, the computation will break down with negative 
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temperature appearing inside the shock wave before a steady structure of the shock wave 
is formed, which is possibly caused by the non-hyperbolic nature of the moment system. 

In this paper, we present a new regularization method which is able to produce smooth 
profile for large Mach number shock waves and low order moment systems. The idea 
originates from the order-of-magnitude method [15] 113]. where the order of magnitude for 
each moment with respect to the Knudsen number is investigated in order to obtain a 
transport system with a specified order of accuracy. Additionally, from the computational 
perspective, a conservative form of moment equations is preferred, so we put this idea into 
the framework of [I] and derive a uniform expression of the regularization terms for all 
moment systems. 

As the first step, we derive the analytical form of the moment equations using the 
same set of moments as in [3J. Once the moment equations are given explicitly, we find 
that only conservative variables and the moments within four successive orders appear in 
each equation. With the help of Maxwellian iteration [8], the order of magnitude can be 
obtained for each moment, and this skill has been used in [10] . The closure of the moment 
system is achieved using a similar skill as in [10] 114]. We approximate the (M + l)-st order 
moments by removing all terms with higher orders of magnitude than the leading order 
in the corresponding equation to get the closed system of all moments with orders lower 
than M. Eventually, a parabolic system is explicitly obtained. 

The resulting regularization term is somewhat complicated and is simplified using 
the technique of linearization for the sake of convenient numerical implementation. As 
in [16], the fluid is considered to be in the vicinity of velocity- free equilibrium states, 
thus the derivatives are small. Dropping the terms which are nonlinear in small values, 
the remaining linear part turns to be very compact. In the ID case, it is obvious that 
the regularization introduces additional diffusion to the M-th order moments. With the 
simplified regularization term, the numerical investigation of the shock tube problem shows 
the convergence to the Boltzmann-BGK equation in moments. And it is numerically 
demonstrated that smooth shock profiles can be obtained for large Mach numbers. 

The layout of this paper is as follows: in Section 2, an overview of Boltzmann-BGK 
model and our discretization of the distribution function is introduced as some prelimi- 
naries. In section 3, the details of the new regularization method are presented. In section 
4, we present two numerical examples to make comparisons between results for different 
moment equations, different Knudsen numbers and different Mach numbers. At last, some 
concluding remarks will be given in Section 5. 

2 Preliminaries 

2.1 The Boltzmann-BGK model 

In the mesoscopic view, the gas can be characterized by the distribution function 
f(t, x, £), where t, x and £ stand for the time, the spatial position and the particle velocity 
respectively, and x,£ £ M. D , D ^ 3. The macroscopic quantities including the density p, 
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the velocity u and the temperature T can be related with / by 

p(t,x)= f f(t,x,£)d£, 
Jr d 

p(t,x)u(t,x) = f Sf(t,x,e)d£, (2.1) 

Jr d 

p(t, x)\u(t, x)\ 2 + Dp(t, x)RT(t, x) = [ \$\ 2 f(t, x, £) d|, 

JR D 

where R is the gas constant. As usual, we use 9(t,x) = RT(t,x) to simplify the notation. 
Since R is a constant, 6 can also be considered as the temperature in the non-dimensional 
case. 

The Boltzmann-BGK model is a simplification of the classical Boltzmann equation, 
which uses a relaxation term instead of the binary collision operator. The Boltzmann- 
BGK equation reads 

^ + i-V x f = \UM~f), (2-2) 
where r is the relaxation time, and Jm is the local Maxwellian defined as 

Mt > x >® = [2*o{t,x)]°/> exp { — m^r ) • (2 - 3) 

By multiplying the Boltzmann equation by (1,£, \£\ 2 /2) T , and integrating both sides over 
~R D with respect to we get the conservation laws as 

k=l 

dm dp A da ik 
^F + ^ + g^ = °' (2 - 5) 

D d9 A dq k A A <9u; 

fc=i «=i j=i j 

where 4j is the material derivative, pij, Uij and q k are the pressure tensor, the stress tensor 
and the heat flux, respectively. The precise definitions are 

dip dtp . . 

-JT = -^7 + u -v x ip, ip = p,Ui,d, 

dt at 



Pij 
Ik 



1 (€i-iH)(€j-Uj)fd£, (t,j pij pd, r p = p6, (2.7) 

n \€-u\ 2 (.£,k-u k )fd£, i,j,k = I,--- ,D, 

2 Jr d 
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where we have used the ideal gas law in p = p9. 

2.2 Discretization of the distribution function 

Suppose x and t are fixed, we expand the distribution function into Hermite functions 
as in |6j 

fit) = E f«HoM, (2-8) 

qGN d 
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where a = (a±, ■ ■ ■ , od) is a .D-dimensional multi- index, and 



v = 2.9 



The basis functions %e,a are chosen as 

D 



1 + 1 / 2 \ 

= J] ^=^^^e ad (^)exp (-^J , (2.10) 
where -f/e Qd is the Hermite polynomial defined by 

He n (x) = (-!)» exp (^) -^exp . (2.11) 

i7e n is assumed to be zero if 7i is negative Thus 1~Lq q is zero if any component of a is 
negative. Some useful properties of the Hermite polynomials are listed in Appendix [A3 It 
has been derived in [3] from these properties that the following relations hold 

D 

fo = P, fe t = 0, ^/ 2ed = 0, i = l,--,D. (2.12) 
d=l 

The stress tensor and heat flux can also be expressed in a simple form: 



Vij = fei+ej , <7jj = 2/ 2ej , 1, j = 1, • • • , D, i^j, 

D 

Qk = 2/ 3efc + J2f^a+e k , k = 1, - ■ ■ ,D. 
d=i 



(2.13) 



In fact, (|2.8|) defines a set of moments M = {fa} a ^D, which will result in an "infinite 
moment system" if we put (j2.8j) into ()2.2|) . In order to get a system with finite number 
of equations, we choose a positive integer M ^ 3 and consider only a subset of Ai which 
contains f a with \a\ ^ M. If we simply force the remaining moments to be zero, the 
Grad-type system associated with the moment set {f a }\ a \^M is obtained. 

In [3], the same set of moments has been used and we have already constructed a 
numerical algorithm for solving the associated Grad-type systems. In the remaining part 
of this paper, we will focus on the regularization of such moment systems. 



3 Regularization of the moment method 

Using a similar technique as in |16j . one possible regularization for the Grad-type sys- 
tems with moments {f a }\ a \^M has been introduced in [3], where a numerical regularization 
algorithm is proposed without deriving the analytical expressions of the regularization 
terms. However, such regularization can cause breakdown of the computation due to the 
appearance of negative temperature while solving the shock structure with Mach number 
Mo ^ 3. This is possibly caused by the non-hyperbolic nature of the moment system. 
In this section, a new regularization method is proposed following the idea of Struchtrup 

H3I- 
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3.1 Maxwellian iteration 



The Maxwellian iteration was introduced by Ikenberry and Truesdell in [8] as a tech- 
nique to derive NSF and Burnett equations from the moment equations. Later in |10j . it is 
used as a tool to analyse the order of magnitude of each moment and to derive equations 
for extended thermodynamics, which is known as the COET method. Below we apply 
Maxwellian iteration to the moment set M, and give a uniform description on the orders 
of magnitude for the moments. 

In order to perform the Maxwellian iteration, we first need to derive the explicit expres- 
sions of the infinite moment equations. This can be done by substituting the distribution 
function / in fj2.2j) with its expansion (|2.8p . After some calculation, both sides of (|2.2p can 
be expanded into Hermite series. Then, we match the coefficient of each basis function 
and then the moment equations can be obtained. The details can be found in Appendix [Bl 
Here, we write the moment equations (|B.8j) as the following form with only one moment 
on the left of each equality: 



, _ _ / df a ,\^dud f 

J"- T \ dt + l^ dt * 



d=l 



D 



+ E 

D 



df a 



d Xj 



2 dt 

d=X 

+ n i^- + («i + 1 ) 



dfa. 



d Xj 



du d 



l£r { d fa-e d - ej + Ujf a -e d + fa + l)f a - ed+ej ) 
d=l i 
D 

2e d -ej + Ujfa-2e d + («j + l)/ a -2e d +e j ) 



d=l 



(3.1) 



lal > 2. 



Note that the cases of \a\ =0 and \a\ = 1 are not included, since when \a\ =0, the 
collision term is zero and the form with /o on the left hand side does not exist, and f a = 
when | a | = 1 following (|2.12p . 

The equation (|3.ip can be taken as an iterative scheme 



(n+l) 



-rG a ( & ] | E N D 



Va E N D and \a\ ^ 2, 



n 



0,1,2,--- 



(3.2) 



with the initial value to be the Maxwellian 



/«S 0) = ^ /r = o, vh^i. 



f (0) 



(3.3) 



During the iteration, /q and f e . are never changed since they never appear on the left 
hand side of (|3.ip . u and 6 do not change with n either. Thus, the operator Q a in ()3.2p 
can be considered as a linear operator according to its analytical form (|3.1|) . For a simpler 
notation, we define the following vectors: 



,(«) 

M 



(fa t)\a\=M, 



= (F 



(ft) 
' ■ 



,(n) 
1 > ■ 



,(«) 
2 J 



(3.4) 



Here is an infinite dimensional vector, but we will reveal that only finite number of 
its components are nonzero. Thus, ()3.2p can be written as 



F (n+1) = _ T g( F (n)^ 



(3.5) 
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and precisely, it is a system as 



fi n+1) = -rG a F^F^F^F^M^) , \a\ > 2. (3.6) 



|a|-3' |a|-2' |o|— 1' |a| ' |a| + l / ' 

Now we start the iteration, and the first two steps will be concretely presented as 
below. 

The first step of iteration As the first step, we start from the initial values and put 
(|3,3p into (|3.2p . Noting that most terms in the right hand side of (|3.2p are zero, we have 

,(D f 1 99 Buj 1 „ \ 

4] = -(^ + ^&J+2'"' V ^J' 



\ dxj dxi 

(1) i 00 (3-7) 



f(l) 

/« = 0, \a\ > 4, 



/^j e , +efc = 0, ^J7^, 



where all /o's are replaced by the density p. Taking r as a small quantity, all moments 
produced in the first iteration are not larger than 0(r). Thus we have 

F W =F (0)_ tF ^\ f (1) =0(1). (3.8) 

The second step of iteration Due to the excessive complexity of the expressions, the 
detailed formulas in the second step of the iteration are not presented while the orders of 
magnitude can be observed. Since Q is a linear operator, we have 

F m = - t q{f^) = -rg(F^) + t 2 g(f {1) ) = fW + r 2 g(F {1) ). (3.9) 

Thus only second order terms are added in this step of iteration. Due to (j3.T|) . the moments 

(2) 

fa with I a I = 2, 3 are not larger than O(r), and the moments with |a| ^ 4 are no larger 
than 0(t 2 ). Furthermore, the moments with \a\ ^ 7 are zeros, which can be revealed by 
(|3.6p and the last equation in (|3.7p . 

Let us go one step closer. For \a\ =3, since 

d df m 

^ ) = --E*-i^ + ---> H = 3 > ( 3 - 10 ) 

j = l 3 

fa^ with | a | = 3 is known to be no less than 0(r 2 ). More precisely, we have 

All = 0(r), / 2 ( 2 +e .=0(r), ff +ej+ek = 0(r 2 ), i + i + k. (3.11) 
For \a\ =4, we have 

D D 



3=1 d=l 



/ W = -^E^^ + - = °(A H=4. (3.12) 
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For \a\ = 5, 6, since D ^ 3, fa^ can be estimated as 



j=l 3 d=l 



(3.13) 



The final conclusion The technique used here can be applied to further iterations 
recursively. It is then found by induction that for any positive integer n, one has 



jrH _ F (n-i) + T n g(F ), p 



0(1). 



(3.14) 



And for \a\ ^ 1 + 3n, /q"^ is zero. This implies that for any a and n, is never larger 
than 0(7~r! a '/ 3 l). Moreover, a careful investigation similar as (|3. 12[) and (|3. 13[) gives 

/(») = O^^H/ 3 !), Va e N D and |a| ^ 4, n > |a|/3. (3.15) 

For |a| ^3, detailed results have been given in (j2. 12|) . (|3.7p and (|3.1ip . 

Remark 1. The equation (|3.14p indicates that for any moment f a , the leading order term 
is never changed once it is obtained. Thus the leading order terms of the stress tensor err- 
and the heat flux qj, can be derived by (|3.7p as 



D 



q k = 2f 3ek + fe k+ 2e d = -^rpO^- 4 
d=l 



fc = l,... ,D, 



®ij — fei+ej 



V: 



2c , 



+ 0(r 2 ), 



J = 1, 



(3.16) 

(3.17) 
(3.18) 



As expected, the Fourier law is deduced as (|3. 16[) . Using (|2.6p . f|3. 18j) can be further 
simplified as 



j j 



dxj 



4 0(r 2 



(3.19) 



where we have used the fact that = 0(t), qk = 0(t). The equations (|3.17p and (J3.19D 
can be written uniformly as 



a 



-2t P 8 



du(i 



dx 



3) 



4 0(r 2 ), i,j = l, 



,D, 



(3.20) 



which is the Navier-Stokes law. It is obvious that the above equations and (|3.16|) yield a 
Prandtl number Pr = 1, which agrees with the common knowledge that the BGK equation 
produces the incorrect Prandtl number 1. This is a validation of the Maxwellian iteration. 

Actually, the Maxwellian iteration can be considered equivalent to the Chapman- 
Enskog expansion, which also gives successive order of the distribution function (see e.g. 
|15|). Here the Maxwellian iteration is much easier to use than the Chapman-Enskog 
expansion, though the latter one is able to give the same results. 
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Remark 2. It is possible in our calculation that some lower order terms will cancel each 
other and only high order terms remain in the iteration. That means some f a with a ^ 4 
may have a smaller order of magnitude o(r^ a l/ 3 T). However, the analysis has depicted 
the essential trend of the variation of the moments' magnitudes, which can be validated 
by comparing (|3. 15[) with the tables in [10]. Meanwhile, the conciseness of (|3.15p is very 
helpful to our later use. 



3.2 The moment closure 



In order to close the moment system, as stated in Section [221 we first choose a positive 
integer M ^ 3 and discard all equations containing the term df a /dt with \a\ > M. Then, 
since f a with \a\ = M + l remains in the system, it is to be substituted by some expression 
consists of lower order moments only. 

This can be done by using (|3.ip and removing the high order terms. With the help of 
(I3.15p . the equation (|3.ip can be reformulated as 



K \d=l d=l / 



D 



dfa 



D 



d=l 3 



+E 

j"=i 
1 00 

+ 2~dx - E ( e f a - 2e <i- e 3 + U jfa-2e d + {Otj + l)f a -2e d + ej ) 
3 d=l 



(3.21) 



D 



+ h.o.t. 



where "h.o.t." stands for high order terms, and it will not be explicitly written later on. 
Note that 



dud _ dud sr^ dud d0 _ 80 sr^ 00 



Putting them into (|3.2ip . we get 

fa=-T 



Edud IdQ . 



D 



3=1 



D 



+ 



9fa 



1 80 



D 



a-2e d 



+E 



— t 

d=l J 



la-e d -ej 



20x~- E^ ] fa-2e d - ej + (aj + l)f a -2e d + e] ) 



d=l 



(3.22) 



(3.23) 



Substituting the material derivatives by the equations (|2.5j) and (|2.6p . (|3.23p is reformu- 
lated as 



L r d=l 7 = 1 J ' 7 = 1 \ J d=l 



dud 

dXj 



D 



d=l 



a-2e d 



D 

E 



df 



D 



d Xj 



+ E 



( dud 



J \ dx, 
l v 



Of, 



a-e d -e. 



1 d0 



(3.24) 
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Recall that a%j and qj have the order of magnitude 0(r), and pij = p5ij + a^. The parts 
containing a^j and qj can also be discarded from (|3.24p . After that, one obtains 




-2e rf 



D 

£ 

j'=i 



9zj 



{Qfa-2e d -ej + (aij + l)/a-2e 



(3.25) 



This equation can be further simplified by coupling it with (|3. 16|) and (|3.20p . and then 
dropping small terms. The final result is 



+ 



dp 

D D 

i=l d=l 



/a 



D 



3/a 



9j(#/a-2e d - ei + («j + l)/a-2e d + ej ) 



(3.26) 



The equation f|3.26j) will be used for all |a| = M + 1, and we ultimately obtain a closed 
parabolic system for the moment set {f a }\ a \^M- 

Remark 3. As is well known, the most serious deficiency of the BGK collision operator is 
that it predicts an incorrect Prandtl number. Therefore, some other models such as the 
ES-BGK [7! and Shakhov [12 j models are proposed as a remedy. Until now, these models 
are known to be very accurate in most cases. All these models have a unified form of the 
collision term: 

Q(f) = u(G-f), (3.27) 

where v is the average collision frequency, and G is some pseudo-equilibrium. The dis- 
cretization of such collision operator has been discussed in [3]. Here we emphasize that 
models with such form can always be easily written as an iteration scheme like (|3.21|) due 
to the existence of the term —/in the collision operator. Thus we can still use Maxwellian 
iteration to analyse the order of magnitude for each moment. 



3.3 Linearization of the regularization terms 

Once the regularization term (|3.26|) is constructed, the system is closed. However, 
recalling that such moment systems are mainly used for computation, it is clear that 
(|3.26p is not concise enough for implementation in numerical simulation. Therefore, we 
are going to linearize (|3.26p and make its expression neater. The similar way has been 
used in [T7] for simplified numerical schemes. 

The linearization will be taken in the neighbourhood of a velocity-free Maxwellian. 
Suppose the radius e of the neighbourhood is small and 

p = po(l + ep), u=\j%eu, 9 = O (1 + e§), 

x = Lex, t = —7=er, f a = PoV tf a for \a\ ^ 1, 
V y o 
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where po and 9q are constants, L is the characteristic length, and variables with : are 
dimensionless of O(l). Thus Oij and qj can be expressed as 



0y — Povoecij, 



Q 3/2 * 

9j = /W e 1r 



(3.29) 



Now we substitute (|3.28p and (|3.29p for the corresponding terms in ()3.26p . After elimi- 
nating the constant factors on both sides, the result reads 



f a = f 



D 



1 y.d(l + ep)(l + e 



D 



9io 



0/o 



J'=l 



<9x 9 - 



D D 



< ~ — i j 1 



-e<j-e,- 



(3.30) 



+ 



i=l d=l 
D D 

EE 



e 2 qj 



After collecting the terms of O(e), (|3.30p is reformulated as 

df, 



(1 + e(9)/ Q _ 2ed -e 3 + {(Xj + l)f a -2e d - 



D 



fa 



-Ed + 



<9f i 



0(6) 



(3.31) 



and the 0(e) term is then simply dropped. Note that one O(e) term is intentionally kept 
in (|3.3ip for the convenience of variable restoration, which is performed by 



fa = pA al/2 ef a = - P0 e l al/2 ef £(1 + 



M/2, 



D 



a/a 



a(poi 



i(|a|-l)/2« 



=er • e (l + rfl) g ^-^ = -rt? £ 



(3.32) 



dxj 



Obviously, (|3.32p is much neater than (|3.26p . and this linearized regularization term is 
used in our numerical examples. 

Remark 4. In the ID case, the regularization term (|3.32p becomes 



fa 



j dfa—ei 

dx 



\a\ = M+ 1. 



And this term is only used in the following term in (|B.8|) 

. a/a+ei 



(ai + 1)- 



<9.r 



lal = M. 



The equations (|3.33p and ()3.34p yields a diffusion on the M-th order term 



d_ 

dx 



(ai + 1)t6 



3x 



|a| = M, 



(3.33) 



(3.34) 



(3.35) 



which reveals the effect of regularization on the Grad-type systems. 
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Remark 5. In the case that M is not a multiple of 3, one may find that the linearized 
regularization term (|3.32j) becomes 0(r^ a ^ 3 ^ +1 ) while it ought to be 0(r^ a l/ 3 T). This is 
caused by using different conceptions of "magnitude" between regularization and lineariza- 
tion. Despite of this, the linear regularization (|3.32p indeed convert the moment system 
to a parabolic one. As known, Grad's moment equations restrict the distribution function 
in a pseudo-equilibrium manifold (see |15|). but the regularization (|3.32p relieves such re- 
striction and allow a small perturbation around the manifold. Although the perturbation 
may not be large enough, it introduces additional flexibility for the moment system to 
agree with the real physics. Also, in our numerical experiments, only slight difference can 
be found between (ET26|) and (|3T32l . 

Remark 6. One may argue that the large moment system is aimed at non-equilibrium 
fluids, and the assumption that the fluid is around a velocity- free equilibrium may lead to 
remarkable deviations. Actually, if we define 

5(0= E ( 3 - 36 ) 

|a|<M-3 

and then linearize (|3.26|) around instead of the Maxwellian, it can be found that the 
linearized result is exactly as (|3.32p . This explains why there is no significant difference 
between the linear and nonlinear regularizations in numerical results when M is large. 

3.4 Comparison with earlier approaches 

In [3], an approach to solve moment system of arbitrary order has been proposed. 
There we use the asymptotic expansion 

f = k + eh + e 2 f 2 + ■■■ (3.37) 

to derive an approximation of /i, where /o is the M-th order Hermite expansion: 

/o = V f a n 9 , a (^) ■ (3.38) 

And the result is 

fi = -r (j£ + ti-V x fX (3.39) 

The similarity between the method above and the current approach is clear. Actually, 
(I3.2ip can be obtained by taking moments on both sides of (|3.39|) and collecting some 
"high order terms". Below we explain the differences between these two methods, which 
are the our major motivation to write this paper. 

The first difference comes from a defect in the theory of the earlier method. In (|3.37p . 
the part f — fo is scaled by a small pseudo-timescale e. However, it is not clear why 
this part can be considered as "small". This is now clarified by the order-of-magnitude 
approach since the magnitude of each moment has been made clear. 

As has been reported in Section [U there exist some circumstances when the earlier 
method fails due to the possible loss of hyperbolicity while the new method does not. 
According to the common theory of the Grad-type methods, this only happens when the 
solution is relatively far away from Maxwellian, which means the "high order terms" that 
have been thrown away in this new approach cannot be simply neglected. However, it is 
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extremely complicated how these terms affect the hyperbolicity of the equations, which 
we have no idea to make clear so far. Only in our numerical experiments, we find that 
dropping these terms alleviates the problem of hyperbolicity. One may argue that the 
new method decreases the accuracy of the approximation, while in our opinion, the loss 
of accuracy can be compensated for by increasing the number of moments. 

Moreover, technical differences originate in designing numerical schemes for these two 
different approaches. In [3], we used (|3.39|) directly as the regularization term. In the 
discretization of (|3.39|) , a direct temporal difference was used to approximate the temporal 
derivative. This is hard to be extended to the second order scheme. And for the part of 
spatial derivatives, it is known now that only three orders of moments contribute to f\, 
which is exactly what we have done in the new method. While in the earlier numerical 
scheme, the difference of the whole distribution function was used for approximation, so 
that all the moments have contribution to this term. Let us demonstrate this point in 
detail in the ID case: the old scheme approximates F = £{V x fo ° n the j-th grid as 



F > ~ 2A^ ' (3 - 40) 



Now, consider the (M + l)-st order moment of Fy. 
Fj,a = C 9j , a [ Ke^aiv^Fjit) exp(K 2 |/2) dVj 

JR D 



a 



(3.41) 

/ -He.^vj) [Uoj+m ~ £i/oj-i(£)] exp(|^ 2 |/2) d Vj 

JR D j 



2Ax 

where Vj = (£ — Uj)/y/9j, \a\ = M + 1, and 



(2ir) D / 2 9[ al+D 

C 9j , a = — r 1 ■ (3.42) 



a l 



Since (uj,6j) / (uj-i,6j-i) and (uj,0j) / Qj+i) in general, the above calculation 

requires projections. Thus all moments of foj+i and foj-i contributes to Fj >a . In the 
new method, we first write the (M + l)-st order moments of F as 

F a = C e , a / Ue^ivj) exp(|^| 2 /2) • i\V x f Q 6.Vj 

JR D 

Q 9f o,a-ei . ^d u d, Q t , f \ 

~ — dx — ^ ^• /o ' Q -^- ei + u iJ0,<*-e d ) (3.43) 

d=l 

1 86 ^ 

,a— 2ed~ ei + Ulfo,a~2e d + («1 + l)/o ,a— 2e d +ei ) j 

X d=l 

and then this equation is adopted to design the numerical scheme. In this expression, it 
is clear that only three orders of moments have contribution to F a . This is the major 
difference between the two methods in the numerical fold, which might lead to large 
deviations in calculating numerical fluxes. The underlying reason of this difference lies 
in different understandings of the regularization term, although such disagreement can be 
eliminated by the refinement of the computational mesh. 

Additionally, the new scheme is more efficient than the earlier one. Currently the most 
expensive part in the algorithm is the projection, which requires to solve an ordinary dif- 
ferential system. The above-mentioned differencing of distributions requires twice of such 
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projection, which is no longer needed in the new framework. Considering the construction 
of a first order numerical scheme, it is found that the computational time can be almost 
halved by such improvement. 



4 Numerical examples 

In this section, two numerical examples of our method for the regularized moment 
equations are presented. In both tests, the global Knudsen number is denoted as Kn, 
and The CFL number is always 0.95. We use the POSIX multi-threaded technique in our 
simulation, and at most 8 CPU cores are used. 

4.1 Shock tube test 

To demonstrate that the new method is applicable to the examples in [3], the first 
example is a repetition of the shock tube problem in [HE]. 
As in [H [T7] , the initial conditions are 



The computational domain is [—1,1.5] and the stop time is t = 0.3. The relaxation 
time t (see (|2.2p ) is chosen as Kn/ p(t,u). The numerical scheme is an improved version 
[5] of the method used in [1], with the regularization term substituted by (|3.32p . The 
improved numerical scheme significantly reduces the computational cost by a large time 
step method and high spatial resolution. Since the BGK model fails to predict the correct 
Prandtl number, we plot the temperature instead of the heat flux below. 

To validate the method with the new regularization term, we compare the results 
produced by the new method with the results in [U [T7] for both small and big Knudsen 
numbers. For small Knudsen number Kn = 0.02, the plot of density for the distribution 
functions when M = 3 and t = 0.3 is presented in Figure [TJ The density profile agrees 
with the result in [H [T7] perfectly. 

For Knudsen number as great as Kn = 0.5, both linearized and non-linearized results 
from M = 4 to M = 15 are computed, which are plotted in Figure [2j Meanwhile, we solve 
the Boltzmann-BGK equation directly using Mieussens' discrete velocity method [5] on a 
very fine mesh grid. One can find that the differences between linearized and the original 
non-linear results are only observable when M is small, which verifies the comments in 
Remark [6) Furthermore, the computational results still converge to the BGK solution 
gradually for both density and temperature, although the convergence rate becomes much 
slower, as can be seen in Figure [2j 

In [18] , it is pointed out that the Grad-type moment system is not globally hyperbolic 
even for 13-moment system. In the case of a large ratio of the density and the pressure, the 
solution leaves the region of hyperbolicity, which leads to strong oscillations and finally a 
breakdown of the computation [18] . Though the ratio of the density and the pressure is as 
large as 7.0 in Figure [21 our method still works well and produces converging results. In 
order to verify this point, an even larger Knudsen number Kn = 5 is investigated. Results 
with M = 3, 6, 9, 12, 15, 18 are considered, and the density and temperature profiles are 
listed in Figure EJ Until M = 18, the regularzied moment method has given a satisfying 
agreement with the discrete velocity model. 




(4.1) 
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Figure 1: Results for the shock tube test with M = 3 and Kn = 0.02. The thick curve 
with the left y-axis is the plot of density. The thin curve with the right y-axis is the plot 
of temperature. 200 grids are used in computation. 



4.2 Shock Structure Problem 

In this section, we carry out the simulation of a steady shock structure with large Mach 
number. The shock structure can be obtained by solving a ID Riemann problem based 
on the Rankine-Hugoniot condition. The left state is 

p r = l, u r = J-M Q , p r = l. (4.2) 



and the right state is 



4M? /5Af n 2 + 3 5M?-1 fl . 

» = wh' — V »— (4 - 3) 

Both states are in equilibrium. After a sufficiently long time, a steady shock wave with 
fully developed structure can be obtained. This example is aimed at the validation of our 
algorithm in high Mach number. For high order moment systems, the current scheme still 
suffers the problem of hyperbolicity. However, the R20 equations (M = 3) turn out to 
be very robust in our numerical experiments. Here we simulate the R20 equations with 
Mq = 1.53, 1.76, 2.05, 2.31, 3.38, 3.8, 6.5, 9.0 and compare the results with the experimental 
data in [2]- The relaxation time is chosen as 

V IhKn 

(4-4) 



2(5-2w)(7-2w) p 

which is the result of the VHS model (see e.g. [3]). The constant u is chosen as 0.72 
as suggested in [2]. The Knudsen number Kn = 1.0 and spatial grid size Ax = 0.1 are 
used in this example. The results for the discrete velocity model are also presented as a 
reference. All the plots are shown in Figure HI The density has been normalized to the 
interval [0, 1]. 

As stated in |llj . in Grad's moment theory, a continuous shock structure exists only up 
to the largest characteristic velocity. For the case of M = 3, i.e. the 20-moment equations, 
Weiss [19] has found that no continuous shock is possible when the Mach number is larger 
than 1.808. In Figure HI the moment system with the new regularization produces stable 
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(a) M = 4 



(b) M = 5 




(e) M = 8 (f) M = 9 

Figure 2: Results for the shock tube test with i^n = 0.5. The dashed lines are the results 
for regularized moment equations without linearization, and the solid thin lines are the 
results with linearization. The dashdot lines are the results of discrete velocity model. The 
black lines denote the density and the gray lines denote the temperature (to be continued). 

and smooth shock structures for much greater Mach numbers. For Mq < 3, the R20 
equations give good agreement with the experimental data, while for larger Mach numbers, 
the profiles are generally correct, although the predicted density is somewhat lower than 
the physical case in the high density region. 
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(g) M = 10 



(h) M = 11 




(k) M = 14 (1) M = 15 

Figure 2: Results for the shock tube test with Kn = 0.5. The dashed lines are the results 
for regularized moment equations without linearization, and the solid thin lines are the 
results with linearization. The dashdot lines are the results of discrete velocity model. 
The black lines denote the density and the gray lines denote the temperature. 

5 Concluding remarks and discussions 

A numerical regularized moment method has been presented. In order to construct the 
regularization term, we first use Maxwellian iteration to determine the order of magnitude 
for each moment, and then approximate the high order moments by eliminating terms 
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(a) M = 3 



(b) M = 6 




(e) M = 15 (f) M = 18 

Figure 3: Numerical results for shock tube problem with Kn = 5. The dashed lines are 
the results for regularized moment equations with linearization, and the solid thin lines 
are the results of the discrete velocity model. The black lines denote the density and the 
gray lines denote the temperature. 
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(a) M = 1.55 



(b) M = 1.76 




(c) M = 2.05 (d) M = 2.31 

Figure 4: The shock structure for different Mach numbers. All lines are density profiles. 
The solid lines are the experimental data, and the dashdot lines are R20 results. The thin 
dashed lines are the results of the discrete velocity model (to be continued). 

with small magnitude. Finally, the approximation is greatly simplified by the strategy 
of linearization. Compared with the regularization in [3], this method not only makes 
it possible to solve high Mach number flow, but also keeps the convergence to the BGK 
solution in moment number. Currently, it is still a challenge to get physical shock profiles 
by this method, and the work on a comprehensive analysis on the shock structure of large 
moment systems is in progress. 
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(e) M = 3.38 



(f) M = 3.8 




(g) M = 6.5 (h) Mo = 9.0 

Figure 4: The shock structure for different Mach numbers. All lines are density profiles. 
The solid lines are the experimental data, and the dashdot lines are R20 results. The thin 
dashed lines are the results of the discrete velocity model. 

Appendix 

A Some properties of Hermite polynomials 

The Hermite polynomials defined in (|2,lip are a set of orthogonal polynomials over 
the domain (— oo, +oo). Their properties can be found in many mathematical handbooks 
such as p]. Some useful ones are listed below: 

1. Orthogonality: / He m (x)He n (x) exp(— x 2 /2) dx = m\V27r5 m ^ n ; 

Jr 

2. Recursion relation: He n+ \{x) = xHe n {x) — nHe n _i(x)\ 

3. Differential relation: He' n (x) = nHe n -\{x). 

And the following equality can be derived from the last two relations: 

[He n (x) exp(-x 2 /2)] / = -He n+1 (x) exp(-x 2 /2). (A.l) 



19 



B Derivation of the moment equations 

In order to derive the analytical form of the moment equations, we need to put the 
expanded distribution (|2.8p into the Boltzmann-BGK equation (|2.2p . The subsequent cal- 
culation will involve the temporal and spatial differentiation of the basis function Tig a (v), 
which will be first calculated as 

D 

D 1 0V 

2 



d=l 



D \a\+D 

(27r)"T0 — 



D 



dvj 

ds 



D 



I He ad+S . d (v d ) exp 



|o| +Dd9„, , N 



(B.l) 



26 ds 



where s stands for t or Xj, j = 1, 2, 3. The partial derivative dv d /ds can be expanded as 

dv d d f£ d -u d \ 1 du d v d d6 



ds 2d ds ' 



(B.2) 



ds ds 

Noting that the recursion of the Hermite polynomials gives 
we conclude from (iBTl) . (|BT2|> and (lB~3l that 

d=l d=l 

Replacing s with i in the above equation, one can have the following expansion of the time 
derivative term in the Boltzmann-BGK equation (|2.2p by some simple calculation: 



(B.3) 



(B.4) 



df 



oeN 1 



dfa. 



& - E (7)f^ + /a 



( 9/ 



E ( It + £ d/ /a ~ ed + 2 at S /a - 2 ^ ) %e ' a - 

QG N D \ <2=1 d=l / 



(B.5) 



Now we consider the convection term. Substituting xj for s in ()B.4p . and making use of 
(|B.3p again, one has 



dx 



3 



df a 



dxj 



dxj 



dxj 



(B.6) 



+ E { e f»-e d ~e 3 + Ujf a -e d + (ttj + l)/a-e d+ej ) 



1 <9# 



+ 2 ^ (Pfa-^d-ej + Ujfa-2e d + (<*j + l)/ Q _2e d+ej ) 



d=l 
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Using /m = foHo,o(v), the relaxation term can be simply expanded as 



-<Jm -/) = -- E fcfre, a {y). 

T T * rf 



(B.7) 



|a|>l 



Finally, we combine (|B.5|) , (|B.6P and (|B.7P and find the ultimate moment equations as 



d=l d=l j 



0f a 



3foL , I . -I \ 



dfa- 



dxj dxj 



d Xj 



D 

£ 

X) T^T { e fa-e d - &j + Ujf a -e d + («j + l)/ a -e d +ej 



(B. 



i de 



D 



+ oThT. ^2 (Qfa-Zej-ej + Ujf a -2e d + («j + l)/a-2e d +e i ) 



where a S N D and 



1 - ^Oa 



1, a = 0, 
0, otherwise. 



(B.9) 
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